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A CLASS OF MONOTONIC SHOCK-CAPTURING DIFFERENCE SCHEMES 
0. Introduction 

A broad class of gas-dynamic problems can be treated within the 

framework of ideal gas theory, i.e., assuming the gas is inviscid and 

thermally nonconductive. The differential equations of gas dynamics 

follow from the laws of conservation of mass, momentum, and energy 

in integral form assuming continuous differentiability of the fluid 

variables, and constitute a system of quasi-linear hyperbolic equations. 

It is well known that in general these may posess discontinuous solutions 

even when smooth initial data are prescribed^. Physically, the presence 

of discontinuities in the solution typically signals the appearance of 

a shock wave. There are two possibilities for correctly describing 

discontinuous flows within the framework of ideal gas theory. The 

first consists of breaking up all the regions In the problem Into 

subregions of smooth flow, which are described by the differential 

equations of gas dynamics, while the discontinuities (boundaries of the 

subregions), are described by conservation conditions. It is important 

to distinguish weak discontinuities (discontinuities in derivatives of 
the fluid variables), tangential discontinuties, and finally, strong 

discontinuities (shock waves). The second approach consists in utilizing 

the conservation conditions in integral form, which allows discontinuous 
solutions. For historical reasons (the differential equations were 

derived earlier), these are referred to as "generalized," iu contrast 
to the continuous classical equations. We note that consideration of 

generalized solutions extends the class of possible solutions, and it 
is therefore necessary to make use of additional considerations in order 

to determine the correct solution. For example the requirement that 
entropy not decrease at a discontinuity permits us to exclude 

Manuscript submitted May 12, 1981. 


1 







rarefaction shocks in the flow of a perfect gas (from a formal point 
of view such a discontinuity is unstable). 

In numerical modelling of gas motions in the presence of shock waves 

both of the above approaches are possible. In the former, the physical 

jump conditions on the discontinuities dividing the regions of continuous 

flow are used both in order to obtain the conditions connecting the fluid 

variables on the two sides of the discontinuity and to determine tha 

motion of the computational mesh points following the discontinuity. 

Changes in the topology of the lines of discontinuity can be followed 

either using a moving curvilinear mesh or by means of an algorithm based 

on a fixed mesh. The first method is most completely worked out in the 

2 

paper by S. K. Godunov et al. while the second has been developed by 

. , 3,4 

Moretti et al. 

The shortcoming of both approaches is the inhomogeneity of the 
resulting difference schemes and consequent complicated structure of lh> j 
numerical algorithm. Both approaches, as a rule, demand an a p riori know¬ 
ledge of the flow pattern during the calculation. 

When it is impossible to know in advance the flow pattern f or when 
it is changing qualitatively with time, it is more convenient to make use 
of "shock capturing" difference schemes, amounting to difference approxi¬ 
mations to the integral form of the conservation laws for every computa¬ 
tional cell. The form of the difference equations does not depend on 
the character of the flow or the position of the possible shocks, end 
therefore such schemes are homogeneous.[Apparently homogeneous 11 
means universal, i.e., not problem-dependent.—DLBJ In this approach 

4 * 

shocks appear as regions of abrupt variation of the fluid quantities. 
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Shock-capturing schemes can also be based on the differential equations. 
For this purpose small corrections are introduced in the equations, 
typically of nonlinear form, analogous to physical viscosity. The 
equations assume parabolic form and permit a smooth solution which 
tends to the solution of the original system as the f, artificial ,r vis¬ 
cosity vanishes. However, difference analogs of the conservation laws 
are satisfied in this technique, generally speaking, only approximately. 
In difference schemes based on the integral form, the mass, momentum, 
and energy conservation laws are satS-^'ed for every computational cell 
to roundoff. Such schemes are describable in terms of fluxes across 
boundaries of the computational cells, and so the conservation laws 
for each computational region are algebraic consequences of the con¬ 
servation laws for cells, i.e., the schemes are conservative,^ 

We note that the first approach, as a rule, distinguishes only 
simple shocks, while the other types of discontinuity are calculated 
as in the shock-capturing method. A difference scheme employed for 
solving practical problems is required to be accurate, that is, it 
must be able to describe flows on comparatively coarse meashes; it must 
be economical; and it should be simple to employ. The accuracy of the 
scheme for calculating smooth flows is determined by the order of the 
approximation. In the neighborhood of a shock, the change in the fluid 
variables is comparable with their magnitude, so that the concept of 
the order of the approximation becomes meaningless. It is well known 
that schemes of the first type lead to smearing of the shock because 
of the strong numerical viscosity. Schemes of second and higher order 
give rise to significantly less smearing of shock runs, but are typically 










nonmonotonic. One of the principle reasons for the appearance of un¬ 
physical oscillations in the neighborhood of a shock is the nonvanishing 
dispersion of the difference scheme, which leads to misrepresentations 
of the form and speed of propagation of physical disturbances, especially 
at short wavelengths. Nonmonotonicity may also be produced by the non¬ 
linearity inherent in real problems.^ 

Artificial viscosity is introduced to suppress unphysical oscil¬ 
lations. Most methods for achieving this, however, do not yield a 
monotonic scheme, and consequently make the localization of shocks worse 
and strongly smear large gradients. This problem becomes especially 
severe in calculations involving strongly shocked flow arising from 
nonsteady interaction of shock waves with one another and with obstacles. 

Recently a number of methods using nonlinear limiters of various sorts 

7-14 

to filter unphysical oscillations have appeared. Many of these are 

highly effective, so that the assertion that r, good schemes of first and 

2 

of higher order smear shock fronts in practically the same degree* 1 
cannot be accepted as correct. 

It is convenient to choose an actual difference scheme in several 
stages. First the dispersion and dissipation properties are studied 
in the linear approximation using the method of differential approxi¬ 
mations or harmonic analysis for the simplest transport equations. It 
is essential to test the schemes selected through linear analysis on 
simple one-dimensional problems having an exact solution. The final 
stage of development of any difference scheme must be its verification 
on model two-dimensional problems. 









In the present paper a nonlinear conservative smoothing procedure 
is proposed in order to achieve monotonicity in explicit schemes of 
higher order. Although this smoothing procedure is applicable to 
arbitrary schemes, the underlying transport algorithm is taken here to 
be second-order MacCormack differencing.^ The results of the linear 
analysis of the dispersion and dissipation properties of this and a 
number of other widely used schemes are presented in Ref. 16. Below, 
the second stage in the choice of an optimal scheme is considered in 
detail. This is conveniently divided into two parts. The first is an 
investigation of the solution of the linear advection equation, which 
facilitates an understanding of the nonlinear properties of the scheme. 

The second part is a calculation of the flow resulting from the evolution 
of a one-dimensional gas dynamic shock. In conclusion, as an illus¬ 
tration of the possibilities of the proposed numerical method, a 
three-dimensional problem involving the interaction of a plane shock 
wave with an obstacle is discussed. 
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1. Construction of second order monotonic difference schemes 


An effective device for constructing monotonic difference schemes 

having second-order accuracy for smooth flows was presented by Boris 
11 12 

and Book. ’ It consists of two stages. In the first a large numeri¬ 
cal diffusion is introduced into the solution, which guarantees the 
monotonicity of the scheme. In the second stage this diffusion Is 
canceled wherever doing so does not introduce new unphysical extrema or 

accentuate existing ones. By means of this approach a family 
of FCT (Flux-Corrected Transport) methods was constructed. 

We consider the simple example of the Cauchy problem for the 
linear advection equation 


where 


3f 

3t 


+ v 


3f 

3x 


0 


( 1 . 1 ) 


f(0,x) » f Q (x). 

The exact solution has the form f(t,x) = f^(x ~ vt). Let be the 
value of the function at the i** 1 grid point of a uniform mesh on the nth 
time level; let 1 + T be the operator which propagates the solution 
from the nth to the (n+l)st level. The form of T depends on the 
particular difference scheme used. We define a diffusion operator D by 


Df i " “i+1/2 " ^1-1/2 " Q6f i+l/2 " Q6f i-l/2 

where 

6f i+l /2 * f i+l “ f r 


( 1 . 2 ) 








Then the first (diffusion) stage can be represented in the fora 

* (1+T+D)f°. 

Now we cancel the diffusion in such a way that in the resulting solu- 
tion no new extrema In comparison with appear, and those already 
present are not amplified. For this purpose in the antidiffusion 
operator A, given by 

M i = " (<? i+l/2 " q? i-l/2 ) ’ 

we limit the fluxes <P i+i y 2 = ® &i ±+l/2 accordin 8 t0 tlie formula 

*1+1/2 " S * ““ t0 ’ 6? i-l/2‘ l^i+l/zl* 

S ^ i+3/2 J>. (1.3) 

where S = sign cp +1/2* 

The transition from the nth to the (n+l)st time level has the fora 

f” +1 = (1+A)?\ = (1+A)(l+T+D)fj . (1.4) 

This technique—the explicit cancellation of the diffusion—leads to 
retention of some of the diffusion in smooth regions [i.e., where the 
limiter (1.3) does not operate and A =? - D], even when T = 0. Conse- 

4 * 

quently Boris and Book proposed two other algorithms: 
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,, phoenIcal ,, 

fj +1 - [(l+A)(1+T) + D]f” (1.5) 

and 

implicit 

f“ +1 - (l+Ar X (l+T-H))f“, (1.6) 

which permit complete cancellation of the diffusion when T = 0 in 

e 

smooth regions. We note that using the algorithm (1*4)-(1.6) leads to 
substantial improvement in the results in comparison with widely used 
schemes. 

In the SHASTA^ scheme the coefficient Q *» 1/8. Subsequently 

Boris and Book investigated in detail the influence of the value of 

magnitude of Q on the dissipative properties of this scheme and deter- 

12 

mined the optimum dependence of Q on the Courant number. However, 
practically speaking, this optimization is meaningful only for a linear 
equation with constant coefficients. For the equations of gas-dynamics, 
as calculations reveal, good results are obtained for Q between 
1/10 and 1/6. 

One can point out three factors responsible for the success of the 
method of Boris and Book. First, the diffusion stage, as shown by 
linear analysis, substantially improves the dispersion properties of 
the scheme. Second, the diffusion is introduced in a conservative 
manner, and, finally, it is done so as to permit excellent localization 
in the neighborhood of the shocks. Evidently it is precisely these 
last two properties which are most important. In Ref. 13 it wds . 
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proposed to interpret the method of Boris and Book as a nonlinear con* 
servative smoothing procedure and to introduce diffusion at the (n+l)st 
time level according to 


f* - (1+A)(l«)(l+T)fj 

or 

f” +1 - (l+A+D)(1+T)f° 


(1.7) 


( 1 . 8 ) 


This procedure fails to improve the dispersive properties and falls 
somewhat short in quality of the results of the algorithms (1.4)-(1.6). 
However, in contrast to the original method of Boris and Book, which is 
essentially one-dimensional in character and which in two- 
dimensional calculations can only be applied via the method of coordi¬ 
nate timestep-splitting of the the difference equations, smoothing in 
the form of Eqs. (1.7)-(1.8) is easily incorporated in essentially 
two-dimensional schemes. For this purpose, after each timestep successive 
smoothing operations are carried out with respect to each of the coordi¬ 
nates. Such smoothing has turned out to be effective in solving a 
whole series of problems.^ ^ 

In the present paper a simpler smoothing procedure is proposed, 
having a local character while retaining the conservative property of 
the difference scheme. A constant diffusion is introduced only in 
regions which are not monotonic, i.e., where the solution has a numerical 
ripple. As a result of this, just as in the method of Boris and Book, 
there is e certain amount of smearing of physical extrema, but for the 
most part only oscillations produced by the numerical scheme are 
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removed. The conservative property is assured by introducing the 
diffusion in the form of fluxes across the cell boundaries. Using 
the notation introduced above, we can write the smoothing at the nth 
time level as 


f“ +1 - (1+T)f" + D*f”, 


(1.9) 


and on the (n+l)st time level in the form 


f^ 1 - (l+D*)(l+T)f^ . 


In these expressions 


D f i " 9 i+i/2 ~ ^i-l/2* 


where 


( 1 . 10 ) 


^1+1/2 


9 i+l/2’ 6f i+l/2 6f i+3/2 K 0 OT 6f i+l/2 6f i-l/2 0 

( 1 . 11 ) 


0 otherwise- 


In order to code this prescription it is convenient to search the entire 
array £^, recording the boundaries across which the flux is nonvanishing f 
and then in a second pass calculate the nonzero fluxes* 

This procedure for introducing artificial viscosity runs substan¬ 
tially faster than the method of Boris and Book, while yielding com¬ 
parable results. Note that, in contrast with the algorithms (1.4)-(1.7), 
smoothing according to the prescription (If8)-(1.10) in regions where 
the fluid quantities vary monotonically in general leaves the solution 
unchanged. 
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2. Solution of the linear transport equation 

In this and the following section the proposed numerical method 

is compared with a number of different schemes, using a one dimensional 

22 

model problem having an exact solution. The schemes of S. T. Godunov , 

15 23 

MacCormack , V. V. Eremin and Yu. M. Lipnitskii , of first, second 

and third order, respectively, are considered and also the non-linear 

12 

scheme of Boris and Book, SHASTX , embodying the algorithm (1.5). These 
schemes are widely employed and have performed well in practice; therefore 
the present analysis permits one to carry out an objective evaluation of 
the proposed method. 

Let us approximate equation (1.1) on a uniform mesh with Ax * 1. 

The computational mesh contains a total of 100 cells. On the boundaries 
periodicity conditions are imposed. The initial conditions are given in 
the form 


fo.5. 


x > L 
1 x :£ L 


where L - 20 and 


cp(x) = 2 


( 2 . 1 ) 


or 

q>(x) - 0.5 + 0.075x. (2.2) 

In Fig.. 1 are shown the results of the calculations after 800 

time steps with Courant number v« 0.2, i.e. for a total running 

time t - Vt/L, after which the initial square wave (2.1) of width 20 . 

cells has been carried across 160 cells. The curves (al), (a2)» (a3), 

correspond to the schemes of Godunov (which coincides in the linear 

case with one-sided differencing), MacCormack (equivalent to the scheme 

«* 

of Lax and Wendroff)t and Ereciin-Lipnitskii; (b). to the SHASTX scheme; (c). 
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(d), (e) correspond to the monotonic MacCormack scheme with smoothing 
according to algorithms (1.7), (1.10), and (1.9). The best results 
were obtained with SHASTX. The SHASTA** algorithm and the monotonic 
MacCormack scheme with the smoothing (1.9) were very slightly inferior 
to this. The smoothing procedures (1.7) and (1.10) give worse results, 
which are comparable, however, with the results obtained from the third-* 
order scheme. The first-order and second-order schemes yielded 
unsatisfactory approximations to the solution. 

Since the solutions of the linear advection,equation using the 

♦ 

algorithms (1.7) and (1.10) were inferior, however slightly, to the algorithm 
(1.9), in the remainder of this section we will present the results of 
calculations carried out using the latter. We must point out, however, 
that these other algorithms will be analyzed below, inasmuch as they 
possess definite advantages in solving real gas-dynamic problems. The 
results of the present section are applicable in problems where linear 
equations are solved, for example, in certain meteorological problems. 

Figure 2 illustrates the "flexibility" of the various schemes in the 
sense of Ref. 2, i.e., the capability of solving problems in which the 
Courant number varies greatly over the computational region. Here 
profiles are shown at t * 3, obtained for Courant number v* 0.6 
(curve 1) and v= 0.2 (curve 2) using first-order (a), second-order (b), 
third-order (c) and the monotonic MacCormack (d) schemesc The 
superiority of the latter is indubitable. The calculation was not 
carried out with SHASTX, which requires that condition v< 0.5 be 
satisfied**. * 

As we remarked in Section 1, the smoothing which removes numerical 


oscillations smears out physical extrema. Consequently, the test 
problem using initial conditions in the form of a triangular profile 
(2.2) constitutes a stringent test of a scheme. Figure 3 shows 
profiles at time t *» 1 (curves 1) and t * 20 (curves 2) for the same 
schemes as in Fig. 2, if v- 0.2. The first-order scheme practically 
"forgets" the initial conditions; the second-order scheme qualitatively 
changes the shape of the pulse as time elapses, while the thid-order 
scheme preserves the "height" of the dlstrubance well. The monotonic 
MacCormack scheme is only negligibly inferior to SHASTX (the dashed lines 





3. The evolution of a gas-dynamic shock 

The problem of the evolution of a gas dynamic shock involves all the 
characteristic features of one-dimenrional ideal gas flow: creation and 
propagation of a shock wave, contact discontinuity and rarefaction fan. 

By comparing the results of the calculation with the analytic solution» 
we can check both the accuracy with which the characteristic flow regions 
are transported and the speed with which they develop*. A comparison 
of the various schemes mentioned above will be presented, based on this 


\ 


problem. 

As was stated earlier, unphysicaloscillations in the neighborhood of 
shocks and steep gradients appear because of the nonvanishing dispersion 


associated with a difference scheme, the nonlinearity of the gas-dynamic 

24 


equations, and also for several other reasons . Besides this, it is 
possible that the various nonlinear smoothing operators may interact 

with one another nonselfconsistently. 

The system of equations has the following form: 


C 1 \ 2 f ' 2 r 

i v x+ l H 


dt = 0, 


f * 2 i t2 r t2 2 \ 2 

I pu dx + I (pu + p)| 
\ *1 ‘ Xl 


dt » 0. 


C 2 \ 2 C 2 

J e I dx + J (e + p )u I 
x h t, x. 


dt * 0. 


'2 , 


Here e = p/(y-1) + pu /2, where y is the adiabatic index. We consider 

<• 

the evolution of a shock with initial pressure ratio p^/p^ m 2 and 
density ratio pj/p^ “ and y « 1.4. In Fig. 4 are shown density 
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profiles obtained using the Godunov (a), Eremin-Lipnitskii (b), 

MacCormack (c) and SRASTX (d) schemes. The exact solution is 
indicated by the heavy continuous line and the numerical results 
up to t « 25 by a fine continuous line. The first-order scheme, as 
one might expect, badly smears the shock front and contact discontinuity 
and represents the rarefaction fan poorly. The Eremin-Lipnitslcii and 
MacCormack schemes are nonmonotonic and yield oscillations in the 
neighborhood of the shocks. SHASTX, which represents the shock front 
very well, gives rise to small undershoots in the density in the neighborhood 
of the contact discontinuity and rarefaction fan. For this case the 
SHASTA scheme, as described by Boris and Book, gives practically the 
same results as does SHASTX. All of the schemes accurately represent 
the pressure profile between the shock wave and the rarefaction fan. 

In Fig. 5 are shown results of calculations using the monotonlc 
MacCormack scheme with various forms of smoothing. Fig. 5a corresponds 
to algorithm (1.7). One shortcoming of this approach lies in the presence 
of undershoots in density in the neighborhood of the shocks and rarefaction 
fans. Very similar results are obtained when one applies the smoothing 
(1.9) to all of the equations in system (3.1) (Fig. 5b). 

In the present work it is proposed to carry out self-consistent 
smoothing of the density, momentum, and energy in regions where 
nonmonotonicity in the density profile exists, using the procedures 
(1.9) and (1.10) described above. The justification for this lies in 
the fact that every stable one-dimensional shock has a density 
discontinuity. In particular, this self -cons is tent smoothing 
introduces no changes in the velocity at a contact discontinuity, since 
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the density and momentum are smoothed in a "related" manner. It also 
somwehat reduces the expenditure of machine time, since in order to 
determine nonmone tonicity, one array of numbers instead of three is 
searched [Cf.(l.ll)j. We note that this type of smoothing cannot be 
implimented using algorithms (1.4) - (1.8) because of their two-stage 
character. 

The density profiles obtained using self-consistent smoothing at 

til ' 

the n and (n+1) st levels are shown in Figs. 5c and 5d. These 

results are clearly superior to all those shown previously. It is 

interesting to observe that the difference in the results obtained by 

T 

smoothing according to the algorithms (1.9) and 1.10) is much smaller 
than in the case of linear advection (Fig. 1). This is Indicative of 
the decisive role played by nonlinear effects in numerical solutions 
of the gas-dynamic equations. 

We pause to draw attention to some peculiarities in the solutions 
we have found. Most of the calculations were run with At « 0.5. 

SHASTA, SHASTX, and the MacCoraack scheme with self-consistent 
smoothing, however, yielded a highly oscillatory solution. The 
oscillations are a consequence of the nonselfconsistency of the smoothing 
and the interaction between the resulting errors in the various fluid 
quantities. A separate paper will be. devoted to a detailed analysis 
of the nonlinear properties of the different schemes, the nonmonotonicj 
resulting from various causes, and ways of contending with the unphysical 
oscillations that result. We therefore will not go inot the reasons for 
the appearance of computational ripples in detail, but content ourselves 
with remarking that they vanish as the time step is reduced. The 1 
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solutions shown in Fig. 4c and 5b were obtained with a time step 
t « 0.25. 

We now compare tins respective running times. Taking as the unit of 
time that required for a single step according to the MacCormack scheme, 
the corresponding running times for the other schemes ate as follows: 
monotonic MacCormack with self consistent smoothing, 1.4; with 
nonselfconsistent smoothing according to algorithm (1.9), 1.5; using 
algorithm (1,7), 2; for the Godunov scheme 2 - 2.5, depending on how 
the "large" quantities are calculated; for SHASTA and SHA;;TX, 3; and 
for the Eremin-Lipnitskii scheme, 2.5. 

We next consider an example in which a strong shock wavs is formed. 
The initial parameters are taken to be the same as in the paper of Boris 
and Book^Cp^/p^ « 480, p^/p^ ” 8 0 y «= 5/3, At: « 0.02). The Eremin- 
Lipnitskii and MacCormack schemes in the absence of ai> adnilioaal 
artificial viscosity do not permit calculation of such a Oov Si;cvr/. 
and the monotonic MacCormack scheme with smoothing according to 
algorithm (1.7) give approximately the s^ame results as SHaSTX, We 
therefore make further comparisons of the results obtained, using the 
method of Godunov, SKASTX, and mono tonic MacCormack .scheme with h* f- 
consistent smoothing given by (1.9) , In Figs. 6 and 7 a) c shown density 
and pressure profiles obtained using these schemes at t.Aw 1 ~ 7. (deshrd 
lines and t « 4 (fine continuous lines). A vertical dashed* dotted lint 
is used to show the initial position of the shock. The vjrtuec: and 
shortcomings of the numerical methods under investigaiion are 
clearly shown ia*Table 1, in which the values of density and relative 
errors at time t ® 4 are shown. In the left column of the table is shown 






the number of computational cells. The Godunov scheme strongly 
smears the shock front and contact discontinuity and approximates the 
exact solution poorly In the region of the rarefaction fan. SHASTX 
is somewhat better than the monotonic MacCormack scheme in representing 
the shock wave and rarefaction fan, but somewhat inferior to it in 
representing the contact discontinuity. In addition, SHASTX gives 
incorrect values in the neithborhood of the initial shock. 

The monotonic MacCormack scheme and SHASTX have good dynamical 
properties. They permit relatively fast determination of the true flow 
pattern (Cf. the Godunov scheme). This last property is especially Important 
in solving nonstationary problems. 
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4. Calculation of three dimensional interaction of a shock wave with 
an obstacle. 

The possibilities of the proposed method are illustrated below for 

the example of reflection of a plane shock wave from a obstacle of 

finite width. The initial stage of the interaction precess is 

investigated, in which the reflected shock wave begins to form and the 

pressure loading on the object is a maximum. We remark that this problem 

is of interest in its own right, inasmuch as experimental determination of 

nonstationary three-dimensional flows encounters serious difficulties. 

The problem is solved in Cartesian coordinates. The mesh is 

uniform. Timestep splitting was employed in solving the difference 
25 26 

equations * which permitted a substantial saving in running time. 

The conservation laws in Eulerian coordinates have the form 


J" P J ^ dv + J 2 (j) pv»n 
v C 1 t l s 

/sv| 2 dv */7# (pw + p) 

v t l s 


r | h r t 2 r _♦ _ 

/ e I dv + I (e + p) V • n 

lt 1 t. s 


ds dt * 0, 


n ds dt = 0, 


ds dt « 0 


We denote by L q (At) the operator which performs a step in the direction o: 


n + 1/3 
ijk 


L (At) £“ 
a ijk. 


L is defined in the following manner: 











n 


n 


7 ijk - f l"k - ST K+s- r ijW. 


n+1/3 1 r, n ~ -^-fF -F )] 

f ijk 2 L£ ijk + f ijk Ah QF ijk i-q,j-r,k-s jJ » 
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q»6,r«6,s«6 
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The smoothing is carried out according to the algorithm (1.9), The 
transverse momentum components were smoothed independently. 

The complete transport operator will look as follows: 

L(2At) * h (At) L (At) 1 (At) L (At) L (At)• 
y x z x y 

Written thus, even if the split operators L ft (At) do not commute, the 
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difference scheme has second-order accuracy over a double timestep 
The stability condition for the split equations is 


At £ min (At )> a - x,y,z, 
a a 


where 


At £ min ( 


Ah 




ijk a 

and a is the speed of sound. This condition turns out bo be 
considerably less stringent than for the unsplit operators: 

Ah 


At £ min ( 


All \ 

Ivl + a/ 3 







Since At^ ( x is the direction of propagation of the shock wave) 

is usually much smaller than At and At , it makes sense to carry out 

y z 

the calculation in the y and z directions with a double timestep. Thus 
using 


L Q (2At) * L Q (At) L a (At), a = y,z. 


we finally obtain 


n+2mAt 


ijk 


L (At)[L (At) L <2At) L (At) L (2At)] 
j x z a y 


,m-l 


n 


. L (At)L (2At)L (At)L (At)f . 

x z x y X j k 


Use of this algorithm instead of 


n+2mAt 


» L m (2At)f 


ijk 


ijk 


permits reduction of the computational volume by a factor of 1.5. 

In the test the computationalmesh consisted of 20 X 20 X 15 cells. 
Calculation of a single case up to t s 1.2 - 1.5 required about 2 hours of 
time on the M4030. Here t « tV /h. where V is the speed of the 

reflected shock wave which arose as a result of the interaction of the 
incident wave with the square object,, and h is the step height. 

A typical flow pattern is shown in Fig. 8. At t ® 1.1 the isobars 
are shown in (a) with separation Ap = 10.0,and the lines of constant 
Mach number (b) at intervals AM = 0.2. The Mach number of the incident 
shock wave M - 5.0 and the ratio of width to height of the step was 
Ji/h 23 5.2. The adiabatic index was y » 1.4 and tha pressure in the 


V 
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undisturbed gas was P^ = 1.0. It is clear that the method permits 
excellent resolution of the shock wave and the flow in the neighborhood 
of the corner points. 

In Fig. 9 are shown the isobars in two cross sections x ** const* 
located equidistantly and parallel to the boundary of the step at a 
distance 0.3 h. The ratio i/h « 2.0 in (a), 3.6 in (b), and 5.2 in (c). 

In Fig. 9a are shown the isobars for the plane case (fc/h * «). 

Figure 10 illustrates the influence of the ratio i/h on the bending 
of the shock wave in the xy plane. Curves 1,2,3 ^correspond respectively 
to i/h = 2/0, 3.6, and 5.2 (the bending is related to the bending for the 
plane case). 

The authors thank Yu. P. Golovachev for his unstinting attention to 
this work and useful discussions, and also V.M. Goloviznin, G*L. Stenchikov^ 
and A.P. Favorskii for valuable suggestions. 
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